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^ . Abstract 

^~ J | We find analytical solutions of the Navier-Stokes equation for the flow 

C/3 . of an incompressible, Newtonian fluid through a porous medium in chan- 

nels with a rectangular and a circular cross section. Nield and others 
• <-H , ° 

CZ3 , have developed an analytical solution for the first case. We express it 

^"»' in an equivalent but a more conventional form. In the second case, we 

i-jH ' find an exact solution when then Forchheimer term can be ignored and 

i *—" *, approximate solution when it cannot be ignored. 

t— I \ 

> . 1 Introduction 

The flow of an incompressible, Newtonian fluid through a porous medium of 
porosity e and permeability K is described by the equation [Zaoli2002 
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where p is the density of the fluid, Ui is the ?th component of its velocity, p 
the pressure, v e the effective viscosity, v the kinematic viscosity, A denotes the 
Laplace operator and the total force is 



&: R = _£^„_^i 



Fi = -— Ui - -^{ujutf^Ui + F M . (2) 

It is a sum of the body force Fu and the effect of the porosity of the medium. 
The geometric function Fo in equation ((2|) depends only on the porosity and the 
nature of the solid particle matrix of the porous medium. Einstein summation 
convention is used in equation @ and others in this paper. For a fully developed 
and steady flow, the left hand side of equation ([T]) vanishes. Further, if the flow 



is driven only by a constant pressure gradient and not an external body force, 
then equation ([TJ becomes 

v e Aui - —m 7 =(u J Uj) 1 ^ 2 u i + ed = 0, (3) 
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where pGi — —dp/dxi. In this paper we solve equation ([1]) in two cases - (a) 
flow between parallel solid plates with a porous medium between them and (b) 
Poiseuille flow in a cylindrical channel stuffed with a porous medium. Nield and 
others Nieldl996 have given an analytical solution of the first case. It expresses 
the coordinate y as an elliptic integral in the velocity u. We shall develop an 
equivalent solution, in a more conventional form, in which u is expressed as a 
Weierstrassian elliptic function of y. We give an exact solution in the second 
case when the term quadratic in u, called the Forchheimer term Zaoli2002 , 
is ignored. If the Forchheimer term is included, we develop an approximate 
analytical solution at the centre of the channel and near its walls. 

2 Flow between two infinite parallel plates 

Consider a steady, fully-developed flow of a fluid along the x\ axis between two 
infinite parallel plates a distance H apart. Let the planes be defined by the 
equations X2 — and X2 = H. If the channel is long enough that the end effects 
are negligible along most of the channel's length then we can express the fluid's 
velocity along the x\ axis as u(x 2 )- Equation ([3]), in this case, becomes 

eG = 0, (4) 

where u" denotes the second derivative with respect to x 2 and pG = —dp/dx\. 
Since the channel is bounded by solid walls, the usual no-slip boundary con- 
ditions, u(x 2 — 0) — and u(x 2 = H) = 0, apply. Let U be the magnitude 
of velocity at the centre of the channel. We convert equation @) in a non- 
dimensional form using the relations, 

u* = - (5) 

x 2 + = — (6) 
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G* = 772- (7) 



and the dimensionless quantities 



Re = ™ (8) 

Da = -^ (9) 

J = ^. (10) 



Re, Da and J are the Reynolds number, the Darcy number and the viscosity 
ratio of the flow. Dropping the * subscript for sake of notation clarity, we get, 

" (U,C: e u -^ u 2 = o (11) 
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Introduce a new variable U\{x2), 



. . / F eRe , 1 / e 

Ul(X2) = V ITdT^ + 2^ JFoReW (12) 



so that equation ([TO becomes, 
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Multiplying equation (fP3"]) by u' 1; the first derivative of u\{y) with respect to y, 
and integrating we get, 

(u' 1 f+4b 2 u 1 -^-u 3 1 +k = 0, (16) 

where fco is a constant of integration. Introduce yet another variable 1*2(2/) = 
(a 3 /36)ui(y) to transform equation (jTB)) to 



( u f,) z =4« 3 -— ^-^ (17) 
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This equation is of the form, 

{u' 2 f = 4u 3 - g 2 u 2 - g 3 , (18) 

where 

fb 2 a 3 \ 
,2 = (— ) (19) 

53 = - W - (20) 

Its solution is 1*2(^2) = p(%2 + k\\ ff2, 53) |Whittakcrl927|, where p is the Weier- 
strassian elliptic function and &i is another constant of integration. Therefore, 
the solution of equation (Ql is, 

F ei?e V V 9 / 36 z J 2F ReVDa 

We find the constants fcp and k\ using the boundary conditions. 



2.1 Finding k 

We shall follow the method suggested by Nield and others Nieldl996 to de- 
termine k . It depends on the fact that at the centre of the channel, the ve- 
locity of the fluid is an extremum (actually, a maximum) and hence its first 
derivative with respect to x 2 vanishes. From the relation between u<i{x<i) and 
u(x2) 1 it follows that whenever u'{x<i) = 0, u' 2 {x2) also vanishes. In that 
case, from equation (fT5|) . it follows that 4u| — g 2 u 2 — 53 = 0. If ei, e 2 and 
e3 are the roots of this cubic polynomial then at the centre of the channel, 
(it2(l/2) — ei)(ti2(l/2) — e2)(tt2(l/2) — 63) = 0. Therefore, one of ei, e 2 or e?, 
is definitely equal to u 2 (l/2), the speed at the centre of the channel. Let us, 
without loss of generality, assume that e 2 = U2(l/2). 
Now, equation (TT8|) can be expressed as, 



U 2 = V 4 "2 - 92U2 - 53, (22) 

or 

nl/2 , r du 2 , s 

dx 2 = / (23) 
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where P$ = u 2 (§). Using the relation between u 2 and u and the boundary 
condition u(0) — 0, it is easy to evaluate P$ to be e/(12JDa). This equation 
can still not be solved because both e 2 and 33 are unknowns. They are however 
related by the following properties of roots of cubic equations, 

ei + e 2 + e 3 = (24) 

eie 2 + e 2 e 3 + e 3 ei = — — (25) 

eie 2 e 3 = ^ (26) 

From these equations, we can easily conclude that 33 = (4e 2 — g 2 )e 2 . We can 
now solve equation (|2"3"]l numerically and get e 2 , and hence g$ and fco, for given 
values of e, i?e, Da, G and J . 

2.2 Finding fei 

We used the boundary condition u(0) = and the fact that in a Poiseuille flow, 
u has a maximum at the centre of the channel to find ko- We will now use 
the second boundary condition u(l) = to find k\. The equation u(l) = is 
equivalent to, 

p (l + k 1 ;g 2 ,g 3 )=P =—^=- (27) 

12J L>a 

This equation always has a number of roots equal to the order of the ellip- 
tic function jWhittakerI927 , 2 in this case. Since u 2 obeys the differential 
equation, 

(p'(z; 52 ,(73)) 2 =4p 3 - 52 p- 5 3, (28) 



where z denotes a complex variable, we can write, 

dz = dp (29) 

V 4 P 3 ~92p- 93 

Further, since the Weierstrassian elliptic function has a pole of order 2 at z = 0, 
we can integrate the left hand side of this equation from p = oo to p = P . 
The corresponding limits of integration on the right hand side are and Pq, 
respectively. Therefore, 

f p=p ° dp 

z(p = P ) - z(p = oo) = z(p = P ) = =, (30) 

J P =oo v 4 P - 92P ~ 93 

or 



i + h= I I =- . ° =, (31) 

?=oo V4p 3 - 32 p - 33 «/p=.Po V 4 P 3 - 92 p - 33 



that is, 



'■■' = - ' - / ^^^^ - / 75 5 ^ (32) 
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The second term on the right hand side is the value of zero of Weierstrassian 
elliptic function. We refer to the papers by Eichler and Zagicr Eich lerl982] or 
Duke and Imamoglu [Dukc2008 for the formulae to evaluate it. The third term 
on the right hand side is an elliptic integral that can be evaluated numerically. 

3 Poiseuille flow without Forchheimer term 

Consider the flow of a viscous, Newtonian, incompressible fluid through a cylin- 
drical pipe of radius R, stuffed with a porous medium of porosity e and perme- 
ability K. Let us, for sake of simplicity, ignore the Forchheimer term which is 
quadratic in velocity. We can do so, if we assume a slow flow. Continuing the 
notation of the previous section, the equation of a fully developed flow is now 
given by (see equation ([3])) 

v e Am - —Ui + eG = 0, (33) 

K 

Although this is a linear partial differential equation, the author is not aware of 
prior work giving its solution for the problem at hand. 

Let us use a cylindrical coordinate system (r, 0, z) with the z axis coinciding 
with the axis of the cylindrical channel. A fully developed flow along the cylinder 
is described by the velocity function u(r)e z , e z be a unit vector along the z axis. 
Equation (|33[) then becomes, 

d 2 u 1 du\ ev 

^ + -rTr)-K U + (G = ^ (34) 



If U is the velocity at the centre of the channel, we can introduce non-dimensional 
variables using the relations, 

u* = £ (35) 

r * = ~R (36) 

G * = w (37) 

and the dimensionless quantities 

Re = — (38) 

Da = § (39) 

J = ^. (40) 

v 

Re, Da and J are the Reynolds number, the Darcy number and the viscosity 
ratio of the flow. Dropping the * subscript for sake of notation clarity, we get, 

u' 

u" H h c - alu = 0, (41) 

r 

where primes now denote differentiation with respect to r and the constants ao 
and Cq are, 



00 = \jjb-a (42) 

co = ^f- (43) 

We will solve equation (14"TT) with the boundary conditions u(r = R) = and 

that u(r) is finite for all < r < R. Let Ui(r) = Co — a^u(r). Equation (|4ip 

transforms to, 

?/ 

< + — - Oo«i = °. ( 44 ) 
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while the boundary condition becomes Ui(r = R) — cq. Equation (1441) is the 

modified Bessel equation of order whose solution is fco/o(ao r ) + kiKo(aor). 

Here, fco and k\ are constants of integration, while 1$ and K q are modified Bessel 

functions of order 0. Since u is finite at r = 0, k\ = 0, while u\(r = R) = cq 

implies fco = co/Io(aoR). Therefore, the solution of (|4"4")l is, 

Mr) = l^R~) IoM (45) 



The solution of (|4i"|) . therefore, is, 



lojapr) \ 



u(r)=ReGDa{l-^^), (46) 

where we have used ao defined in equation (J42]) to simplify notation. 



3.1 Asymptotic form of solution 

We shall now examine the behaviour of the solution P5|) as the porosity of the 
medium tends to 1, that is in the limiting case of no porous medium. To that 
end, write PrJ|) as, 

. . RcGrDcL ,,,,,, . * 

U M = T7 - m ( J o(aoi?) - Jo(oor)) , (47) 

lo(aoK) 

The series form of the modified Bessel function of order m is (formula 9.6.10 of 
|Abramowitzl964j ). 



k— 

Therefore, 

_ fleGDa // e \ (fl 2 -r 2 ) / e V* (R A - r A ) \ 

u[r) ~ io(a R) \\jm) 4r(2) + V lm) ~¥2W{^ + ■ ■ ■ ) : { ' 

The absence of porous medium is characterised by the limits e — > 1, J — > 1 and 
Da — > co, or in terms of our constants, ao — > 0. 

lim u(r) = ^(i? 2 -r 2 ), (50) 

which is the solution of equation (l4"Tj) in the limits ao — >• and e — > 1. 

4 Poiseuille flow with Forchheimer term 

We shall now consider the Forchheimer term in the equation equation of flow. 
That is, we will attempt to solve the equation, 

fi; F 1 

u e Au u ^=u 2 + eG = 0, (51) 

K y/K 

subject to the boundary condition u(r — R) — 0. After introduction of non- 
dimensional variables as before, we can transform it to, 

u' 
u" H h c - a 2 w - b u 2 = 0, (52) 

where the constants ao and Co are same as before while bo is defined as, 

; F ° Re (**\ 

b a = (53) 

J\l Da 



Equation (l5Tj) is similar to equation (|4Tj) except for the term bou 2 . It is also 
similar to equation (|11|) except for the additional term u'/r. However, unlike 
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equations (|41l) and (fTTj). it is not easy to obtain a solution of this equation. 
We shall therefore attempt an approximation. We shall use the solution of 
the problem without the Forchheimer term as a guide to approximation. In 
particular, we shall examine the behaviour of v! jr and u" at the centre of the 
channel and at its walls. From equation (|46|) . 

•w - ^"('"Is) ,54) 

, (ReGDa)a 

U{r) = I (a R) /l(a ° r) (55) 

u = jRzGDaK f oM _hM\ (56) 

J (a i?) V a o»' / 

These relations follow from, 

^-{ x PI p {x)) = x%_i(a:) (57) 

J_ p (z) = I p (x) (58) 

From the series representation of I m (x), we can easily show that, 

u"(r) 
lim - \ . = 1 (59) 

r-rt) {u'{r)/r) 

r^fl (u'(r)/r) V /l(oo-R) / 

Typically, ao is a large number (because Da is usually very small). Therefore 
equation ([50)1 suggests that near the walls of the channel, the second derivative 
of velocity dominates u'(r)/r. However, at the centre of the channel, both terms 
are of equal magnitude. 

4.1 Approximate solution near the channel's wall 

We can ignore the u'/r term near the wall of the channel. Equation (|52|) can 
then be approximated as 

u" + Co — OqM — boil 2 = (61) 

A transformation u\ = u^/b~a + a 2 /{2^/ba) gives, 

Multiplying by u[, and integrating with respect to r gives, 

«) 2 + 2 (c + |0 U! - H„f + fc = 0, (63) 



where fco is a constant of integration. Introducing ui — 6ui gives, 

(u' 2 ) 2 = 4^-12^0 + ^L\u 2 -mk , (64) 

whose solution is a Weierstrassian elliptic function p(r + fci; (72,53), with 

to = 12(co + ^) (65) 

53 = 36/c (66) 

and k\ being a constant of integration. The approximate solution at the centre 
of the channel is, 

« r) _£ p ( P+fciU („+£),,«,)_£ (67) 

The constant k\ can be computed in a manner similar to section[2j However, fco 
will have to be chosen so to adjust the solution near the wall to that mid-way 
between centre and the wall. 

4.2 Approximate solution at the centre of the channel 

We observed that near the centre of the channel, the second derivative of the 
velocity is of the same order of magnitude as the u'(r)/r term. Therefore, we 
assume that the equation of flow is approximated by 

u' 
2— + c - alu - b u 2 = (68) 

r 

Using the transformation U\ = Uy/bo + a 2 /{2\/b^)^ we get 

2ui 1 u o \ 2 



It can be further simplified to, 



c o + 7T - < = (69) 



dUl =^rdr, (70) 



u\-d\ 2 

where the new constant do is defined as, 

Equation (|70|) can be readily integrated to give, 

^ = -wM( d ^ +k °)}-i' (72) 

where kg is a constant of integration. There is no boundary condition available 
at the centre to determine fco and like the identically named constant in the 
previous section, it too will have to be determined by adjusting the solution at 
the centre to that valid mid-way between centre and the wall. 
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